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ABSTRACT 

The oceanic turbulent boundary layer is a critical region to 
understand for oceanic and atmospheric prediction. This thesis 
answers two fundamental questions :( i ) what is the response of the 
ocean mixed layer system to transient forcing at the air sea 
surface? (ii) what is the necessary time and space resolution in an 
ocean mixed layer model to resolve important transient responses? 

Beginning with replication of de Szoeke and Rhines' (1976) 
work, additional physical processes were added to include more 
realistic viscous dissipation and anisotropy in the three- 
dimensional turbulent kinetic energy (TKE) budget. These refinements 
resulted in modification of de Szoeke and Rhines' findings. 
Firstly, TKE unsteadiness is important for a minimum of 10 s seconds. 
Secondly, viscous dissipation should not be approximated as simply 
proportional to shear production. Thirdly, entrainment shear 
production remains significant for a minimum of one pendulum-day. 

The required temporal model resolution is dependent on the 
phenomena to be studied. This study focused on the diurnal, 
synoptic, and annual cycles, which the one-hour time step of the 
Naval Postgraduate School model adequately resolves. The study of 
spatial resolution showed unexpectedly that model skill was 
comparable for 1 m, 10 m and even 20 m vertical grid spacing. 
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I . INTRODUCTION 

A. STATEMENT OF PROBLEM 

The oceanic turbulent boundary layer or mixed layer is a 
critical region to understand when studying the atmosphere or 
ocean. The ocean is the primary source of moisture fluxed 
into the atmosphere, it is a significant source of heat fluxed 
into the atmosphere in many regions, and it is a moderator of 
atmospheric extremes. To accurately model the atmosphere or 
ocean, a realistic representation of the turbulent mixed layer 
is required. This thesis will address two fundamental 
questions: (i) what is the response of the ocean mixed layer 
system to transient forcing at the air sea surface?; (ii) what 
is the necessary time and space resolution in an ocean mixed 
layer model to resolve important transient responses? With 
the answers to these questions, air-ocean models will be a 
step closer to being realistically coupled. The result will 
be better predictions for both atmosphere and ocean. 

B . METHOD 

To accomplish this task, the work of de Szoeke and 
Rhines(1976) on wind driven mixing was replicated and then 
generalized to include surface heating and cooling. Their 
study in "Asymptotic regimes in mixed layer deepening" did not 
include physical processes associated with surface heat flux 
and anisotropy in the three-dimensional turbulent kinetic 
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energy (TKE) budget. With the addition of these critical 
physical processes, a more generalized study of the short time 
scale dependence of the mixed layer is possible. The study 
can be extended past the single day scale of de Szoeke and 
Rhines to the synoptic and annual time scales by shifting the 
method of solution from a Runge-Kutta numerical scheme to a 
FORTRAN-based gridded model. With the question of the minimum 
time scale resolved, new assumptions about mixed layer physics 
can be made which allow testing of spatial resolution. The 
Naval Postgraduate School (NPS) mixed layer model developed by 
Garwood (1977) was modified to test the effect on model skill 
of varying the vertical grid size using observed 
meteorological forcing and changing thermal structure at 
weather station Papa (50° N, 145° W) . 
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II. THEORY, MODEL EQUATIONS, AND HYPOTHETICAL SOLUTIONS 



A. AIR OCEAN COUPLING 

The surface layer of the ocean responds to variations in 
wind stress, heat flux, moisture flux, and several other less 
significant factors. These factors have been examined in 
depth since Ekman's (1905) ground laying work. 

B. REPLICATION OF DE SZOEKE AND RHINES (1976) 

Initial modeling efforts for this study focused on 
reproducing and then building on the efforts of de Szoeke and 
Rhines (1976) as discussed in their paper "Asymptotic regimes 
in mixed-layer deepening." Niiler (1975) proposed the basic 
terms in the turbulent kinetic energy (TKE) equation studied 
by de Szoeke and Rhines: 

( — N 2 h A - u? — — (1 -cos ft) + C n uf h 2 ) =2m n uf h (1 ) 

2 f 2 0 dt 0 

A B CD 

Term A represents buoyant damping of TKE by entrainment, term 
B the shear production of TKE by entrainment, term C the spin 
up or storage of TKE, and term D the near-surface wind shear 
production of TKE minus viscous dissipation. Solar heating 
and surface heat flux were not included. Term C, the energy 
required to spin up the turbulence intensity, was not 
originally included by Niiler. 
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To mathematically define the components of the de Szoeke 
and Rhines model, the friction velocity or u. must first be 
defined. The friction velocity (equation 2) is defined as 
the square root of the surface wind stress divided by the 
density of water. 



u 





P 0 



( 2 ) 



The Brunt-Vaisala or buoyancy frequency (N) is a measure of 
vertical stability. A value of 211/600 sec' 1 was used for 
initial model testing and in the hypothetical cases . The 
mixed layer depth is h. The Coriolis parameter (f) is defined 
as 2Q Z where Q z is the vertical component of planetary 
rotation. A nominal mid-latitude value of lxlO' 4 sec' 1 was 
used for f. In equation (1) C 0 and m 0 are dimensionless tuning 
coefficients for TKE spin up and net wind-shear production, 
respectively . 

Buoyant damping of TKE by entrainment, term A, is the 
rate of decrease in turbulence due to the lifting upward of 
denser water originating below the mixed layer. The mixed 
layer is made more stable by the entrainment of colder denser 
water, with a commensurate increase in potential energy. 
Shear production due to entrainment, term B, is the conversion 
of mean kinetic energy to turbulent kinetic energy due to 
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shear across the base of the mixed layer. This term also 
results in a net gain in turbulence. Turbulence spin up, term 
C, is the TKE build up required before entrainment can occur. 
The wind stirring minus viscous dissipation, term D, 
represents the net input of energy into the mixed layer due 
to wind stress, x 0 , in equation (2) . In de Szoeke and Rhines' 
model, viscous dissipation is also proportional to u. 3 because 
surface heating and cooling are neglected as sources or sinks 
of TKE, and dissipation of TKE produced by entrainment shear 
production is not accounted for. 

As a first step in demonstrating and understanding the 
physics that drive mixed layer deepening, de Szoeke and 
Rhines' results were replicated. Equation (1) was solved for 
the rate of mixed layer deepening, W e = dh/dt . The result is 
referred to as the W e equation, since the vertical change in 
h is due only to entrainment (no advection) . The W e equation 
was formulated as a MATLAB function, and solved using an 
intrinsic Runge-Kutta solution routine, ODE45. The MATLAB 
function is restricted to continuous functions. This 
restricts the model to deepening scenarios . 

De Szoeke and Rhines simplified equation (1) by analyzing 
the time evolution of the mixed layer deepening process and 
determined which terms could be neglected during various time 
regimes. Analytic solutions to the equation were then 
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possible in those specific time regimes. Figure 2.1a is a 
reproduction of de Szoeke and Rhines results, and 2.1b depicts 
the results of a Runge-Kutta solution to equation (1) . 

De Szoeke and Rhines showed that initially shear 
production due to entrainment (term B) can be ignored, that 
turbulence spin up (term C) dominates, and that the mixed 
layer deepens at a rate proportional to u.t. This condition 
holds for the first 100 seconds of deepening. Next, surface 
wind stirring (term D) dominates. Deepening continues at a 
rate proportional to t 1/3 . Term C begins to build in 
importance between 20 and 120 minutes, with the deepening rate 
proportional to t 1/2 . Coriolis (term B) becomes a factor as the 
mean flow direction is turned away from the wind direction. 
After twelve hours, term B is no longer dominant. Deepening 
returns to a wind driven scenario and a rate proportional to 
t 1/3 . 
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Figure 2.1a de Szoeke and Rhines solution. Solid curve: 
ratio D/A; dashed curve: ratio B/A; dotted curve: C/A . 




Figure 2.1b MATLAB ODE45 solution, replicating de Szoeke and 
Rhines . 
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C. IMPROVED MODEL PHYSICS TO INCLUDE EXPLICIT TKE 

In the second model, the basic W e equation is the same as 
in the first model, but three additional equations are 
included: the mixed layer temperature equation and the two 
horizontal mixed layer momentum equations. This allows the 
explicit calculation of the variables of temperature and 
velocity. The full set of equations are: 



dr _ Op 
dt p C h 

p 



A TW 



Su 

dt 



— f v + 



h 

&UW 



(3) 



ph 



dv - X y 
= - fu + - 



h 

A VW 



(4) 



(5) 



dt 



ph 



0 3 k2 a ^ 0 h 3 

2m n u t h - 

3h = PS 

dt (^-W 2 h 4 -u. 4 (1-cos ft) +C 0 u. 2 h 2 ) 



-W(z = -h) 



( 6 ) 



In this more complete model with explicit momentum and heat 
equations, the ability to add steady state upwelling was 
included in the W e equation. This model is identical to the 
original de Szoeke and Rhines model if the upwelling term, 
W(z=-h), is set to zero, and if the net surface heat flux(Q 0 ) 
is neglected. However, when Q 0 is nonzero in the second 
model, Q 0 directly affects both the temperature equation (3) 
and the mixed layer depth equation (6), and it indirectly 
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affects momentum because of the W e terms in equations (4) and 
(5) . In the momentum equations, t >; and t, are components of 
wind stress. AT is the change in temperature across the mixed 
layer interface (Figure 2.2); similarly, AU and AV are the 
changes in the velocity components across the interface. 




Figure 2.2 Hypothetical temperature profile. 

The second model's solution for equation (6) was 
identical to de Szoeke and Rhines when Q c and W(z=-h) were set 
to zero. The three additional equations allowed the explicit 
calculation of mixed layer temperature and velocity, but did 
not change the W e solution. The model was initialized with 
both horizontal velocity components equal to zero and an SST 
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of 17° C with a linear temperature gradient with depth, having 
N = 211/600 seconds. A steady wind stress of 0.1 N/m 2 was 
applied in the i x direction. In 30 hours, the SST decreased 
to 16.3 C (Figure 2.3) as the cooler sub-mixed layer water was 
entrained. The horizontal velocities, u and v, respond to the 
steady state wind stress in the form of forced inertial 
motion. The period is dependent on latitude (2n/f) but a 90 
degree phase shift between u and v is constant. 

The third model is a MATLAB version of the NPS mixed 
layer model; it represents a significant addition of TKE 
physics to the first two models. Equations for diurnal 

heating and cooling and three for energy transformation were 
added, and the W e equation was reformulated. The resultant 
equations are: 

Q 0 = ~Q a + £>n (max (0, eos-^-) ) ( 7 ) 

0 a a 86400 

— - = 2m, u? + U 2 W +2ra, (£-3 u 1 ) Je-— (m.E 2 + m.fhE) (8) 

dt 3 e 2 v 3 ' 1 5 

— - = V 2 W +2®. (E-3 v 1 ) y/E- — (m. £ 2 +m.fhE) (9) 

dt 3 

o w 2 h — 7 , — 2 nn\ 

= -aghA TW + 2 m. ( £ - 3 w 2 ) Je - — (ra, £ 2 +m.fhE) ( 10 ) 

dt 6 pC 2 v 3 1 5 

K p 

w = J ~ 7 1 — — (id 

e aghAT+£ 
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Figure 2.3 Solution for temperature and velocity. 
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( 12 ) 



— =W - W(z = -h ) 
dt e 

A Q 0 (t) equation (7) was designed to have a diurnal structure 
with a net zero heat flux. The energy transformation terms 
(m 2 ) redistribute energy among components, representing 
pressure-strain interactions. The variable E is the summation 
of the squares of the turbulent velocity fluctuations in the 
x, y, and z directions, e=u 2 + v 2 + w 2 • In equation (8), m 3 u, 3 
is the near surface wind shear production term. The summation 
of the three m 1 E 3/2 + n%fhE terms found in equations (8), (9), 
and (10) is the viscous dissipation. The m 3 u, 3 term and the 
summed m 1 E 3/2 + m 5 fhE terms are approximated in the de Szoeke 
and Rhines equation by the m 0 u. 3 h term. The m 2 (E-3u 2 ) E 1/2 terms 
are pressure redistribution of TKE among the components, in 
which the turbulence tends toward isotropy by pressure-strain 
interactions. The aghATW e term in equation (10) is the 
buoyant damping of TKE by entrainment and corresponds to de 
Szoeke and Rhines N 2 h 4 term. The aghQ 0 /pC p term in equation 
(10) is the buoyant damping of TKE by surface heating (Q 0 ) , 
and has no corresponding term in equation (1), although it was 
included in equation (6) of the second model. Equation (11) 
is the entrainment rate as a function of the ambient TKE. 
Equation (12) is the rate of change of the mixed layer due to 
entrainment, W e , and vertical mean advection by W(z=-h). 
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Initial runs of this model were conducted with a constant 



wind stress but no heating or cooling. These solutions were 
similiar to the results of models one and two, but the mixed 
layer depths were significantly shallower (Figure 2.4). This 
is as expected and is attributed to energy being expended for 
turbulence spin up and dissipation in the layer as well as 
entrainment. Two differences between the first and third 
models are significant: (i) calculation of dissipation and 
(ii) role of inertial oscillations. 

Term D, in the de Szoeke and Rhines model, near-surface 
wind shear production of TKE minus viscous dissipation, is an 
over simplification. The viscous dissipation should not be 
represented as directly proportional to u. 3 h. The dissipation 
resulting from the shear production due to entrainment is not 
represented at all. Figure 2.5 is analogous Figure 2.1, and 
it shows the spin up of model terms. The ratio of D3/A3 is 
near-surface wind shear production of TKE minus viscous 
dissipation over buoyant damping of TKE by entrainment. The 
ratio of B3/A3 is shear production of TKE by entrainment over 
buoyant damping of TKE by entrainment. The ratio C3/A3 is the 
turbulence spin-up term over buoyant damping of TKE by 
entrainment. The most obvious difference is in the ratio 
D3/A3. Total dissipation is much larger in the third model 
than with the approximation used in the first model because 
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depth (m) 




Figure 2.4 Comparison of mixed layer depths between model two 
and model three. 
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Figure 2.5 Model three solution. The ratios represent the 
same physical processes as Figures 2.1a and 2.1b. 
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the first model estimated viscous dissipation as proportional 
to the wind stress only. The third model calculates the 
dissipation for all three TKE components based on the 
explicitly calculated TKE in the mixed layer. Between 10 3 
seconds and 4xl0 4 seconds, there is more viscous dissipation 
than wind shear production, because of the addition of the 
dissipation due to shear production. 

The significance of term B decreases as the mixed layer 
deepens. In the first model, the effect of inertial forcing 
is negligible after 6xl0 4 seconds. In the third model, term 
B3 remains significant four times longer than for model one. 
If the dissipation due to shear entrainment is neglected, the 
amplitude of term B3 decreases, and term A3 and D3 dominate 
(Figure 2.6) after lxlO 5 seconds. 

D. MODEL RESOLUTION IN THE TIME DOMAIN 

Models one through three all use a Runge-Kutta solution 
with a variable time step to solve the mixed layer equations . 
The numerical scheme maximizes the time step to solve the set 
of equations without exceeding specified error tolerances. In 
mixed layer physics, dh/dt is not a continuous function when 
mixed layer shallowing occurs. For this reason, the Runge- 
Kutta solution can only be used during periods of mixed layer 
deepening. To study realistic cases, a model that can shallow 
as well as deepen was required. 
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Figure 2.6 Model three solution with viscous dissipation 
Coriolis parameter m 5 =0. 
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The NPS mixed layer model was chosen. The NPS model is 
a FORTRAN-based gridded numerical model with a fixed time step 
of one hour. Using a fixed time step to simultaneously solve 
differential equations introduces errors. To determine the 
adequacy of the one-hour time step, the NPS model was compared 
to the third model using mixed layer depth (Figure 2.7) and 
temperature (Figure 2.8) as the basis of comparison. The NPS 
model does not have a term for shear production of TKE by 
entrainment (term B3) but vertical mixing by dynamic 
instability occurs when the Richardson number (Ri) is less 
than the critical value, Ri < 0.25, following the procedure of 
Adamec et al. (1981) . These physical processes can then be 
selectively included or neglected for comparison purposes. 

The NPS model deepens too quickly during the first time 
step. This is because the NPS model assumes fully spun-up 
steady state turbulence. This is accomplished by solving the 
steady state equations (13), (14), and (15), or by setting 
the left hand side of equations (8), (9), and (10) to zero. 

0=V 2 W +2m 0 (E-3v 2 ) Je- — {m,E 2 +m,fhE) (13) 

e z v ^ i j 



cxghQ 2 4 

0--cxghh.TN + 2m, (£-3 w 2 ) Je- — (m^E 2 +m,fhE) (14) 

e p C 2 3 1 5 

p 

? 2 

0=2m 3 u? + U 2 W e + 2m 2 (E-3u 2 ) JE- (m 1 E 2 +m 5 fhE) (15) 
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Figure 2.7 Mixed layer depth comparison between the NPS mixed 
layer model with Richardson Number deepening turned off and 
model number three with Term B3 turned off. 



19 



Depth (m) 



0 



T 



T 



5- 



NPS 

— Model 3 



10 - 



15 L 

16,1 



16.2 16.3 16.4 16.5 16.6 16.7 

Temperature (degrees C) 



Figure 2.8 Mixed layer temperature comparison between the NPS 
mixed layer model and model number three, at time = 8 hours. 
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This assumption allows instantaneous deepening to occur 
without having to spin up the turbulence. Depending on the 
rate of fluctuation in the forcing, this may not be a valid 
assumption. Figure 2.9 depicts a portion of the C3/A3 ratio 
from Figure 2.5. Examination of the results depicted in 
Figure 2.9 reveals that the unsteady term does not approach 
zero at 10 4 seconds as proposed by de Szoeke and Rhines but at 
10 5 seconds. At 0.36xl0 4 seconds, this assumption results in 
a maximum 14 percent error in TKE generation for any change in 
forcing. If the forcing is constant or slowly changing this 
error is negligible. Figure 2.9 indicates that the unsteady 
term persists and oscillates in excess of 3 days. 

The temperature profiles differ for two reasons. The 
first is related to mixed layer depth. The NPS mixed layer 
model depth is shallower and therefore warmer. Conversely, 
model three has a deeper mixed layer depth and is cooler. The 
second reason is related to model thermal structure. Model 
three has a continuous function for temperature versus depth. 
The NPS model is a gridded model and assumes a step like 
function for temperature. The difference in thermal structure 
between the models is evident in Figure 2.8 in the region 
below the mixed layer. 

To further test the adequacy of the one hour time step, 
the NPS model was run with a 15 minute time step. The results 
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Figure 2.9 Ratio of Model three turbulence spin-up over 
buoyant damping of TKE by entrainment (C3/A3) . 
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are depicted in Figure 2.10. After the second time step of 
the one hour model, the mixed layer depths of the two NPS 
models are parallel in time with an offset of approximately 
0.5m. The one-hour time step of the NPS mixed layer model is 
sufficient for the proposed study. This conclusion was 
reached because the error induced by the relatively large time 
step is small and explainable. After initial spin up, the 15 
minute time step does not result in significant resolution 
improvement. Additionally, the subsequent forcing data to be 
used has a sampling period of three hours. Interpolating to 
time steps less than an hour was therefore judged to be 
inappropriate . 
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Figure 2.10 Mixed layer depth comparison using the one-hour 
time step and the 15 minute time step NPS mixed layer model, 
and model number three. 
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III. REALISTIC SOLUTIONS 



A. IMPORTANCE OF VERTICAL RESOLUTION 

In any geophysical process model, the resolution required 
is dependent on the scale of the physical phenomena to be 
represented. The annual fluctuation in the mixed layer depth 
at station Papa is on the order of 150 m. The daily 
fluctuations vary between order 10 m and order 100 m, 
depending on time of year and the synoptic forcing. Unlike 
other mixed layer models, the NPS mixed layer model used here 
employs a floating grid point assigned to the mixed layer 
depth to increase the model's resolution. This feature 
becomes more significant as grid size increases. 

The choice of vertical grid spacing (Az) is governed by 
several factors. The two factors considered in this thesis 
are model computational speed and resolution. These two 
factors are not independent of each other. Smaller grid 
spacing requires increased computer time, but resolves more 
features; larger grid spacing is faster, but has reduced 
resolution. Model requirements will determine which one of 
these two factors is more important. In the past, embedding 
a mixed layer within an ocean general circulation model or 
coupling a well-resolved upper ocean model to a global 
atmospheric model has been prohibitive because of the number 
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of grid levels required to represent the ocean. The following 
sections address the question of required vertical grid 
spacing and the accuracy of larger grid sizes. 

B . MODEL FORCING 

The atmospheric forcing utilized was calculated from 
observations by a Canadian weather ship in the vicinity of 
50° North latitude 145° West longitude, nominally referred to 
as station Papa. Three hourly weather observations were 
taken, consisting of wind speed and direction, air and sea 
surface temperature, dew point, and fractional cloud cover. 
The three-hourly observations were interpolated to hourly 
values. Aerodynamic bulk formulas derived by Large and Pond 
(1982) and modified by Oberhuber (1993) were used to compute 
the momentum, sensible, and latent heat fluxes, correcting for 
atmospheric stability. Mechanical bathythermographs (MBT) were 
also taken at irregular intervals at station Papa, varying 
from less than an hour to several days. The MBT data was 
manually digitized at 5 m intervals. The MBT data are used 
in this study to initialize the mixed layer model and then for 
subsequent verification. Model runs were conducted for 
various periods between 1961 and 1969. In addition to being 
one of the longest running continuous oceanic time series, the 
data set is ideally suited to testing a one dimensional model 
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because the horizontal currents are weak and advection due to 



upwelling is generally believed insignificant (Tabata, 1965) . 

C. EFFECT OF VERTICAL GRID SIZE 

The 1 m grid NPS mixed layer model has a minimum mixed 
layer depth of 1 m, or one Az . Changing the grid spacing 
raised the issue of what the minimum mixed layer depth or 
ceiling should be. When Az was increased to 10 m and 20 m, 
a 1/4 Az minimum mixed layer depth was found to be optimal. 
The 1/4 Az ceiling resulted in a 30% decrease in mean error in 
temperature over the one-Az ceiling. In the presence of 
strong heating and no wind stress, ceilings smaller than 1/4 
Az resulted in the mixed layer temperature becoming 
unrealistically hot. The same unrealistic temperatures were 
also obtained in the 1 m grid case by raising the ceiling 
above one Az . 

D. VERIFICATION WITH OBSERVATIONS 

1 . Model Predicted SST vs Bucket Temperature 

Sea surface temperature is the only state variable for 
which model output and station Papa observations can be 
continuously compared. To accomplish this task, the model 
mixed layer temperatures were compared to the interpolated 
bucket temperatures. The mean difference and standard 
deviation were used for comparison. 



27 



The mean difference = 



E„ <SST" lc -SST 



bucket • 



n 



The mean difference calculation was used to evaluate the 
effect of varying the minimum allowable mixed layer depth as 
described in section III C. The results presented in Appendix 
A provide the justification for selecting 1/4 Az as the 
minimum ceiling for Az's of 10 m and 20 m. The eight annual- 
period model runs commenced in mid-March, near the time of 
maximum mixed layer depth, and ran for 365 days. The results 
of the mean differences were unexpected and raised questions 
about optimal spatial resolution. In every case but one, the 
error for model runs with 1 m grid spacing was greater than 
for both the 10 m and 20 m Az cases. This implies that the 
coarser grid is more accurate. The standard deviation showed 
the same trend. The model runs depicted in Figure 3.1 are 
characteristic of the other years. Spectra were calculated 
for 1 m, 10 m and 20 m grid resolution ( Figures 3.2, 3.3, & 
3.4, respectively). The annual cycle was not resolved due to 
the spectral computational length of the time series; however, 
the diurnal cycle is evident. The diurnal cycle was best 
resolved by the 10 m grid, then by the 20 m grid, and lastly 
by the 1 m grid. 

In an attempt to understand why better model performance 
is achieved with coarse resolution, model tuning sensitivity 
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Figure 3.1 Observed bucket temperatures and model output for 
1966. 
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C*C*Day 




Figure 3.2 Spectra for 1 m grid spacing for 1966, using 
tunning coefficients ?3 = 4.6 and AM3 = 4.5. 
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C*C*Day 




Figure 3.3 Spectra for 10 m grid spacing for 1966, using 
tunning coefficients P3 = 4.6 and AM3 = 4.5. 
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• C*C*Day 




Figure 3.4 Spectra for 20 m grid spacing for 1966, using 
tunning coefficients P3 = 4.6 and AM 3 = 4.5. 
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tests were conducted. Tuning coefficients P3 and AM 3 in 
equations (13), (14), and (15) were varied, and the model 
rerun for each year. Tuning coefficient P3 equals mj/mj, and 
regulates the effect of the Coriolis in the dissipation term; 
larger P3 results in stronger dissipation. Tuning coefficient 
AM3 equals m 3 and regulates the efficiency of shear 
production. Larger AM3 results in increased shear production. 

Year-long model runs for eight years were made for each 
P3 and AM3 pair. The values ranged for 3 to 8 for P3 and 1 to 
9 for AM3 . The standard deviation for each run was calculated 
and the results averaged. Figure 3.5 depicts the average 
error for the 1 m grid spacing cases; Figures 3.6 and 3.7 
depict the average error for the 10 m and 20 m grid spacing 
respectively. A well-defined region of minimum error is 
evident for the Az equal to 10 m and 20 m cases. This same 
region. is not as clearly defined in the Az equal to 1 m case. 

Optimal tuning parameters were chosen which represent the 
(P3, AM3) region of minimum error. The model was rerun for 
five different time frames and for seven different tuning 
coefficient pairs. The first was a repeat of the annual 
cycle. The next four were 10 day runs, one per season. This 
was done not only to elucidate differences in model accuracy 
among the seasons, but also to determine if the 1 m Az model 
improved over the larger grid models for diurnal and synoptic 
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Figure 3.5 Average standard deviation error (° C) for 1 m 
grid spacing. 
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Figure 3.6 Average standard deviation error ( = C) for 10 m 
grid spacing with a ceiling of 2.5 m. 
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Figure 3.7 Average standard deviation error (° C) for 20 m 
grid spacing with a ceiling of 5 m. 
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scales . 



The standard deviations were calculated and 



tabulated. Appendix B is the complete listing of results by 
tuning factor pairs. 

The results of the in-depth tuning study are both 
clarifying and conclusive. The vertical grid spacing accuracy 
is dependent on time of year and length of model run (see 
table 1) . During the spring and summer when the forcing is 
dominated by the diurnal cycle, the larger Az's are more 
accurate. During the synoptically forced deepening period 
of the fall, the 1 m grid spacing more accurately represents 
the bucket temperatures. In the 48 cases where the grid size 
is a factor when using 10 m vice 1 m Az, an increase in 
standard deviation error of 0.06° C resulted, and an increase 
of 0.05° C resulted in the 35 cases when 20 m Az were used. 
The winter accuracy is generally independent of grid spacing. 

2 . Bucket vs Mechanical Bathythermograph SST 

The NPS mixed layer model uses BT's to initialize the 
model. After initialization, surface meteorological 
observations are used to drive the model . The surface 
observations were taken at 3 hourly intervals, but the MBTs 
were deployed at irregular intervals. For initial conditions, 
the model uses the BT with the closest Julian date to the 
start date of the run. In general the relative vertical 
structure of the surface layer will be correct, but a random 
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bias in the absolute temperature will be induced. 
Additionally, a phase shift is also likely, as the model 
starts at 00 GMT and the BT may or may not be at 00 GMT. 
Figure 3.8 demonstrates the bias between bucket and BT 
temperatures and the irregular sampling interval of the BTs. 
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Table 1. Values represent the number of relative minima in 
standard deviation of the error. The last four columns 
account for model runs were the standard deviations of the 
error are equal. 
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Temperature (degree C) 




Figure 3.8 Example of bucket temperatures vs mechanical 
bathythermograph SST at station Papa. 
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IV. SUMMARY AND RECOMMENDATIONS 



This study has evaluated the response of the ocean mixed 
layer to transient forcing and the time and space resolution 
required to accurately model upper ocean variability. This 
study began with replication of de Szoeke and Rhines (1976) 
work, using a Runge-Kutta MATLAB solution method. Additional 
physics were added to their model including: surface heat 
flux, improved viscous dissipation, and anisotropy in the 
three-dimensional turbulent kinetic energy (TKE) budget. The 
surface heat flux was added for completeness but was not 
extensively tested in evaluating the most rapid transient 
responses with a less than diurnal period. The refinement of 
viscous dissipation and TKE budget resulted in modification of 
de Szoeke and Rhines' findings. First, the unsteady TKE term 
is important for a minimum of 10 5 seconds rather than only 10 4 
seconds. Next, viscous dissipation should not be approximated 
as being proportional simply to u, 3 h. This parameterization 
fails to take into account buoyancy flux due to surface heat 
flux and entrainment. The time scale of the turbulence also 
plays a role in dissipation, and is not previously included. 
Finally, shear production due to entrainment remains 
significant for a minimum of one pendulum-day instead of the 
one-half pendulum-day of the de Szoeke and Rhines' study. 
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Expanding the investigation to model predictions longer 
than a day, the NPS FORTRAN mixed layer model was employed. 
The required temporal resolution of mixed layer models is 
dependent on the phenomena to be studied. This study focused 
on the diurnal, synoptic, and annual cycles, for which the 
one-hour time step of the NPS FORTRAN model is more than 
adequate to resolve. The one-hour time step is well matched 
with the larger vertical grid spacings . Further research is 
required to determine whether a one-hour time step is too 
coarse when used with 1 m grid spacing. If a time step 
shorter than one hour is to be used, the unsteady term may 
need to be added to the NPS FORTRAN model . 

Finally, the question of spatial resolution was 
addressed. The initial hypothesis for this research was that 
smaller grid spacing would yield better results than larger 
grid spacing. The results reported in Appendix A and Table 1 
were therefore unexpected. With the exception of the winter 
time frame, the 10 m and 20 m grid spacing generally equaled 
or out-performed the 1 m grid spacing. Use of the 20 m grid 
spacing during the winter resulted in an rms error of 0.05° C, 
which is smaller than the 0.1° C precision of the bucket 
temperatures. Overall, the NPS FORTRAN model skill was 
comparable for 1 m, 10 m, and even 20 m vertical grid spacing. 
Therefore, the conclusion of this thesis is that use of 10 m 
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and 20 m vertical grid spacing in global ocean models and 
coupled air-ocean models to represent the mixed layer is 
accurate and justified. 
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APPENDIX A. ANNUAL MEAN OF PREDICTED TEMPERATURE ERROR (DEGREES C) FOR 
VARIOUS MINIMUM ALLOWABLE MIXED LAYER DEPTH 



Year 


Az 


1 Az 


0.5 Az 


0.25 Az 


1.01 m 


1961 


1 


-1.5698 


-1.5537 


-1.5257 


-1.5644 




10 


-1.0555 


-1.0425 


-1.0272 


-0.9973 




20 


-1.1308 


-1.0693 


-1.0514 


-0.9948 


1962 


1 


-0.4187 


-0.3665 


-0.3338 


-0.3959 




10 


0.1393 


0.1937 


0.2296 


0.2859 




20 


-0.0428 


0.0795 


0.1181 


0.2257 


1963 


1 


-0.2282 


-0.2014 


-0.1808 


-0.2274 




10 


0.3013 


0.3154 


0.3331 


0.3608 




20 


0.1765 


0.2451 


0.2680 


0.3266 


1964 


1 


-1.2219 


-1.1993 


-1.1838 


-1.2228 




10 


-0.7613 


-0.7451 


-0.7333 


-0.7150 




20 


-0.8379 


-0.8106 


-0.7974 


-0.7618 


1965 


1 


-1.4082 


-1.3979 


-1.3831 


-1.4124 




10 


-0.9455 


-0.9213 


-0.9069 


-0.8899 




20 


-1.0784 


-1.0427 


-1.0276 


-0.9896 


1966 


1 


-0.6312 


-0.6080 


-0.5760 


-0.6246 




10 


-0.1622 


-0.1347 


-0.1143 


-0.0815 




20 


-0.3113 


-0.2362 


-0.2094 


-0.1407 


1967 


1 


-0.5077 


-0.4958 


-0.4786 


-0.5077 




10 


-0.0159 


0.0013 


0.0163 


0.0416 




20 


-0.1604 


-0.0947 


-0.0762 


-0.0289 


1968 


1 


-0.4024 


-0.3775 


-0.3423 


-0.4074 




10 


0.0859 


0.1252 


0.1555 


0.2064 




20 


-0.0529 


0.0047 


0.0373 


0.1382 
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APPENDIX B. RMS TEMPERATURE ERROR (DEGREES C) VERSUS MODEL AZ AND 
TUNNING PARAMETERS P3 AND AM3 





P3=3 . 5 


AM3=2 












Winter 






Spring 








1 


10 


20 


1 


10 


20 


61 


0.27 


0.27 


0.27 


0.28 


0.27 


0.27 


62 


0.17 


0.18 


0.18 


0.21 


0.20 


0.20 


63 


0.20 


0.20 


0.20 


0.17 


0.17 


0.17 


64 


0.14 


0.14 


0.14 


0.11 


0.08 


0.09 


65 


0.13 


0.13 


0.13 


0.22 


0.19 


0.17 


66 


0.14 


0.13 


0.13 


0.10 


0.09 


0.09 


67 


0.35 


0.42 


0.46 


0.11 


0.09 


0.09 


68 


0.19 


0.18 


0.18 


0.07 


0.07 


0.07 




Slimmer 






Fall 








1 


10 


20 


1 


10 


20 


61 


0.63 


0.61 


0.59 


0.21 


0.22 


0.24 


62 


0.74 


0.71 


0.69 


0.44 


0.50 


0.54 


63 


0.65 


0.66 


0.68 


0.24 


0.33 


0.26 


64 


0.66 


0.58 


0.61 


0.33 


0.42 


0.48 


65 


0.72 


0.64 


0.60 


0.45 


0.54 


0.44 


66 


0.72 


0.72 


0.71 


0.08 


0.09 


0.10 


67 


0.72 


0.70 


0.63 


0.16 


0.19 


0.21 


68 


0.65 


0.63 


0.63 


0.17 


0.17 


0.17 




Year 














1 


10 


20 








61 


1.21 


0.49 


0.51 








62 


0.43 


0.95 


0.67 








63 


1.26 


0.53 


0.59 








64 


1.15 


0.48 


0.50 








65 


1.03 


0.42 


0.43 








66 


0.54 


0.58 


0.52 








67 


1.06 


0.54 


0.64 








68 


0.98 


0.45 


0.41 
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P3=4 . 6 


AM3=4 


.5 
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Spring 








1 


10 


20 


1 


10 


20 


61 


0.27 


0.27 


0.27 


0.27 


0.27 


0.27 


62 


0.18 


0.20 


0 .19 


0.21 


0.20 


0.21 


63 


0.14 


0.14 


0.14 


0 . 11 


0.08 


0.09 


65 


0.14 


0.14 


0.14 


0.18 


0.18 


0.17 


66 


0.14 


0.14 


0.14 


0.10 


0.09 


0.09 


67 


0.35 


0.41 


0.47 


0.10 


0.09 


0.09 


68 


0.18 


0.18 


0.18 


0.07 


0.07 


0.07 
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Fall 




* 
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10 


20 


1 


10 


20 


61 


0.63 


0.62 


0.59 


0.20 


0.22 


0.24 


62 


0.74 


0.72 


0.69 


0.42 


0.47 


0.54 


63 


0.67 


0.67 


0.69 


0.22 


0.32 


0.30 


64 


0.67 


0.59 


0.59 


0.26 


0.39 


0.48 


65 


0.72 


0.65 


0.61 


0.38 


0.48 


0.38 


66 


0.73 


0.72 


0.73 


0.09 


0.09 


0.10 


67 


0.70 


0.69 


0.65 


0.16 


0.19 


0.21 


68 


0.65 


0.64 


0.64 


0.17 


0.17 


0.17 




Year 














1 


10 


20 








61 


1.33 


0.87 


0.91 








62 


0.53 


0.43 


0.39 








63 


1.39 


0.92 


0.99 








64 


1.24 


0.82 


0.90 








65 


1.09 


0.66 


0.76 








66 


0.68 


0.35 


0.36 








67 


1.29 


0.93 


1.03 








68 


1.09 


0.66 


0.70 
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P 3 = 7 . 5 


AM 3=5 












Winter 






Spring 








1 


10 


20 


1 


10 


20 


61 


0.27 


0.27 


0.27 


0.27 


0.27 


0.27 


62 


0.20 


0.20 


0.20 


0.21 


0.21 


0.21 


63 


0.20 


0.20 


0.20 


0.18 


0.17 


0.18 


64 


0.14 


0.14 


0.14 


0.10 


0.08 


0.09 


65 


0.14 


0.14 


0.14 


0.18 


0.18 


0.17 


66 


0.14 


0.14 


0.14 


0.10 


0.09 


0.09 


67 


0.47 


0.50 


0.54 


0.11 


0.09 


0.09 


68 


0.18 
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0.07 


0.07 
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20 
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10 


20 


61 
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0.61 


0.58 


0.22 


0.25 


0.25 


62 


0.73 


0.72 


0.68 


0.50 


0.55 


0.59 


63 


0.66 


0.66 


0.68 


0.27 


0.37 


0.27 


64 


0.65 


0.59 


0.60 


0.38 


0.50 


0.54 


65 


0.70 


0.64 


0.60 


0.46 


0.56 


0.43 


66 


0.72 


0.71 


0.71 


0.09 


0.11 


0.11 


67 


0.69 


0.69 


0.62 


0.18 


0.21 


0.22 


68 


0.64 


0.63 


0.63 


0.17 


0.17 


0.17 




Year 














1 


10 


20 








61 


0.85 


0.51 


0.50 








62 


0.54 


0.79 


0.72 








63 


0.89 


0.59 


0.66 








64 


0.77 


0.44 


0.44 








65 


0.65 


0.40 


0.42 








66 


0.49 


0.67 


0.64 








67 


0.89 


0.65 


0.71 








68 


0.75 


0.57 


0.60 
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P 3 = 7 . 5 


AM 3=6 
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0.27 
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0.20 
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0.21 


0.21 


0.21 
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0.20 


0.20 


0.20 


0.18 


0.18 


0.18 


64 


0.14 


0.14 


0.14 


0.10 


0.08 


0.09 


65 


0.14 


0.14 


0.14 


0.18 


0.18 


0.16 


66 


0.14 


0.14 


0.14 


0.10 


0.09 


0.09 


67 


0.48 


0.50 
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0.10 


0.09 


0.09 


68 


0.18 


0.18 


0.18 


0.07 


0.07 


0.07 
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1 


10 


20 
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0.62 


0.61 


0.59 


0.22 


0.25 


0.25 


62 


0.73 


0.72 


0.68 


0.50 


0.55 


0.59 


63 


0.66 


0.66 


0.68 


0.28 


0.37 


0.27 


64 


0.66 


0.59 


0.60 


0.38 


0.50 


0.54 


65 


0.69 


0.64 


0.60 


0.45 


0.55 


0.41 


66 


0.72 


0.71 


0.72 


0.09 


0.10 


0.11 


67 


0.69 


0.69 


0.62 


0.18 


0.21 


0.23 


68 


0.64 


0.64 


0.63 


0.17 


0.17 


0.17 




Year 














1 


10 


20 








61 


0.89 


0.56 


0.55 








62 


0.53 


0.72 


0.68 








63 


0.92 


0.64 


0.71 








64 


0.79 


0.50 


0.54 








65 


0.68 


0.41 


0.44 








66 


0.49 


0.60 


0.59 








67 


0.96 


0.68 


0.78 








68 


0.79 


0.60 


0.62 









50 





P 3=7 


AM 3=7 












Winter 






Spring 








1 


10 


20 


1 


10 


20 


61 


0.27 


0.27 


0.27 


0.27 


0.27 


0.27 


62 


0.20 


0.21 


0.21 


0.21 


0.21 


0.21 


63 


0.20 


0.20 


0.20 


0.18 


0.17 


0.18 


64 


0.14 


0.14 


0.14 


0.10 


0.08 


0.09 


65 


0.14 


0.14 


0.14 


0.17 


0.18 


0.17 


66 


0.14 


0.14 


0.14 


0.10 


0.09 


0.09 


67 


0.44 


0.48 


0.53 


0.10 


0.09 


0.09 


68 


0.18 


0.18 


0.18 


0.07 


0.07 


0.07 




Summer 






Fall 








1 


10 


20 


1 


10 


20 


61 


0.62 


0.61 


0.59 


0.22 


0.24 


0.25 


62 


0.73 


0.72 


0.69 


0.47 


0.53 


0.58 


63 


0.66 


0.66 


0.69 


0.26 


0.37 


0.28 


64 


0.66 


0.59 


0.59 


0.35 


0.48 


0.53 


65 


0.70 


0.65 


0.61 


0.42 


0.52 


0.39 


66 


0.73 


0.72 


0.72 


0.09 


0.10 


0.11 


67 


0.69 


0.69 


0.63 


0.18 


0.21 


0.23 


68 


C . 64 


0.64 


0.63 


0.17 


0.17 


0.17 




Year 














1 


10 


20 








61 


1.01 


0.66 


0.70 








62 


0.49 


0.59 


0.57 








63 


1.01 


0.81 


0.82 








64 


0.90 


0.58 


0.71 








65 


0.78 


0.46 


0.56 








66 


0.51 


0.47 


0.45 








67 


1.05 


0.84 


0.87 








68 


0.85 


0.61 


0.67 









51 





P 3=8 


AM 3=4 












Winter 






Spring 








1 


10 


20 


1 


10 


20 


61 


0.27 


0.27 


0.27 


0.27 


0.27 


0.27 


62 


0.20 


0.20 


0.20 


0.21 


0.21 


0.21 


63 


0.20 


0.20 


0.20 


0.18 


0.17 


0.18 


64 


0.14 


0.14 


0.14 


0.11 


0.08 


0.09 


65 


0.14 


0.14 


0.14 


0.19 


0.18 


0.17 


66 


0.14 


0.14 


0.14 


0.09 


0.09 


0.09 


67 


0.50 
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